perm filename SUBSLM.OLD[DRW,LCS] blob
sn#108366 filedate 1974-12-13 generic text, type T, neo UTF8
00100 SUBROUTINE SS
00200 COMMON X(100),Y(100),N,X1(512),Y1(512),S(100),K
00300 COMMON/SLM/A(2,200),B(2,200),C(2,200),I(200),I1(200),
00350 1 I2(200),KT,TT
00360 DATA L1/1/,L2/2/,L0/0/
00400 C FIND MAX AND MIN POINTS AND SET SLOPES
00500 ICOUNT=0
00600 710 DO 100 K=2,N-1
00610 KP=K-1
00620 KQ=K+1
00700 IF(Y(K).GT.Y(KP).AND.Y(K).GT.Y(KQ)) GO TO 20
00800 IF(Y(K).LT.Y(KP).AND.Y(K).LT.Y(KQ)) GO TO 20
00900 IF(X(K).GT.X(KP).AND.X(K).GT.X(KQ)) GO TO 30
01000 IF(X(K).LT.X(KP).AND.X(K).LT.X(KQ)) GO TO 30
01100 GO TO 100
01200 20 S(K)=0
01300 ICOUNT=ICOUNT+1
01400 I(ICOUNT)=K
01500 GO TO 100
01600 30 S(K)=1001
01700 ICOUNT=ICOUNT+1
01800 I(ICOUNT)=K
01900 100 CONTINUE
02000 ICOUNT=ICOUNT+1
02100 I(ICOUNT)=N
02200 C FIND POINTS SUCH THAT X(K)=X(K+1) OR Y(K)=Y(K+1)
02300 DO 900 K=2,N-1
02400 IF(Y(K).EQ.Y(K+1)) I2(K)=1
02500 IF(X(K).EQ.X(K+1)) I2(K)=1
02600 900 CONTINUE
02700 C SET THE SLOPES BETWEEN
02800 IA=1
02900 IX=0
03000 110 IF(IX.EQ.ICOUNT) GO TO 200
03100 IX=IX+1
03200 IA=IA+1
03300 120 IF(IA.EQ.I(IX)) GO TO 110
03310 KP=IA+1
03400 S(IA)=(Y(KP)-Y(IA-1))/(X(KP)-X(IA-1))
03500 IA=KP
03600 GO TO 120
03700 200 CONTINUE
03800 C SET BEGIN AND END SLOPES
03900 CC SZ=1002001.
04000 K=1
04100 CALL S101
04200 CALL S102(X(1),Y(1),X(2),Y(2),S(1),S(2))
04300 IF(S(K).EQ.-1001.) S(K)=1001.
04400 IF(S(K+1).EQ.-1001.) S(K+1)=1001.
04500 K=N-1
04600 CALL S101
04700 CALL S102(X(N),Y(N),X(N-1),Y(N-1),S(N),S(N-1))
04800 IF(S(K).EQ.-1001.) S(K)=1001.
04900 IF(S(K+1).EQ.-1001.) S(K+1)=1001.
05000 C RESET THE SLOPES FOR STRAIGHT LINES
05100 DO 610 K=1,N-2
05200 U1=SLM2(L0)
05300 U2=SLM2(L1)
05400 IF(ABS(U1-U2).GT..0001) GO TO 610
05500 I1(K)=1
05600 I1(K+1)=1
05700 S(K)=U1
05800 S(K+1)=U1
05900 S(K+2)=U1
06000 610 CONTINUE
06100 C ADD POINTS
06200 K1=N-1
06300 DO 300 K=1,K1
06400 IF(I1(K).EQ.1) GO TO 300
06500 IF(I2(K).EQ.1) GO TO 300
06600 FLG=0.
06700 IF(S(K).EQ.0..AND.S(K+1).EQ.0.) FLG=1.
06800 IF(S(K).EQ.1001..AND.S(K+1).EQ.1001.) FLG=2.
06900 IF(FLG.EQ.1..OR.FLG.EQ.2.) GO TO 205
07000 CC SZ=1002001.
07100 CALL S101
07110 KP=K+1
07200 U=(Y(K)-Y(KP)-S(K)*(X(K)-X(KP)))/(S(KP)-S(K))
07300 AX=X(KP)+U
07400 AY=Y(KP)+S(KP)*U
07500 XA=X(K)
07600 1610 XB=X(KP)
07700 IF(XA.LE.XB) GO TO 202
07800 XA=X(KP)
07900 XB=X(K)
08000 202 YA=Y(K)
08100 YB=Y(KP)
08200 IF(YA.LE.YB) GO TO 204
08300 YA=Y(KP)
08400 YB=Y(K)
08500 204 CONTINUE
08600 IF(AX.GE.XA.AND.AX.LE.XB.AND.
08700 1 AY.GE.YA.AND.AY.LE.YB) GO TO 290
08800 205 K1=K1+1
08900 DO 210 K2=K+1,N
09000 K3=N+K+1-K2
09010 KP=K3+1
09100 X(KP)=X(K3)
09200 Y(KP)=Y(K3)
09300 I1(KP)=I1(K3)
09400 I2(KP)=I2(K3)
09500 210 S(KP)=S(K3)
09600 IF(FLG.EQ.1..OR.FLG.EQ.2.) GO TO 280
09700 N=N+1
09710 KP=K+1
09720 KQ=K+2
09800 Z0=X(K)**2-X(KQ)**2
09900 Z1=Y(K)-Y(KQ)-S(K)*(X(K)-X(KQ))-((X(K)+X(KQ))/2.
10000 1 -X(K))*(S(K)-S(KQ))
10100 Z2=X(K)**3-X(KQ)**3-3*(X(K)-X(KQ))*X(K)**2
10200 1 -1.5*(X(K)+X(KQ))*Z0+3*X(K)*Z0
10300 AZ=Z1/Z2
10400 BZ=(S(K)-S(KQ)-3*(X(K)**2-X(KQ)**2)*AZ)/
10500 1 (2*(X(K)-X(KQ)))
10600 CZ=S(K)-3*AZ*X(K)**2-2*BZ*X(K)
10700 DZ=Y(K)-AZ*X(K)**3-BZ*X(K)**2-CZ*X(K)
10800 X(KP)=-BZ/(3.*AZ)
10900 Y(KP)=AZ*X(KP)**3+BZ*X(KP)**2+CZ*X(KP)+DZ
11000 S(KP)=3*AZ*X(KP)**2+2*BZ*X(KP)+CZ
11100 IF(Y(KP).LT.YB.AND.Y(KP).GT.YA) GO TO 290
11200 X(KP)=(X(KQ)+X(K))/2.
11300 Y(KP)=(Y(KQ)+Y(K))/2.
11400 S(KP)=0.
11500 GO TO 290
11600 280 N=N+1
11700 IF(FLG.EQ.2.) GO TO 282
11800 S(K+1)=2.*(Y(K+2)-Y(K))/(X(K+2)-X(K))
11900 GO TO 284
12000 282 S(K+1)=(Y(K+2)-Y(K))/(2.*(X(K+2)-X(K)))
12100 284 X(K+1)=(X(K+2)+X(K))/2.
12200 Y(K+1)=(Y(K+2)+Y(K))/2.
12300 290 IF(S(K).EQ.-1001.) S(K)=1001.
12400 IF(S(K+1).EQ.-1001.) S(K+1)=1001.
12500 300 CONTINUE
12600 C FIT THE POINTS WITH N-1 CURVES
12700 305 DO 400 K=1,N-1
12750 KP=K+1
12800 IF(I1(K).EQ.0) GO TO 310
12900 A(1,K)=0.
13000 A(2,K)=0.
13100 B(1,K)=X(KP)-X(K)
13200 B(2,K)=Y(KP)-Y(K)
13300 C(1,K)=X(K)
13400 C(2,K)=Y(K)
13500 GO TO 400
13600 310 CALL S101
13700 B(1,K)=2*(Y(KP)-Y(K)-S(KP)*(X(KP)-X(K)))/
13800 1 (S(K)-S(KP))
13900 B(2,K)=S(K)*B(1,K)
14000 A(1,K)=X(KP)-X(K)-B(1,K)
14100 A(2,K)=Y(KP)-Y(K)-S(K)*B(1,K)
14200 C(1,K)=X(K)
14300 C(2,K)=Y(K)
14400 400 CONTINUE
14500 C CALCULATE 512 POINTS
14600 M=N-1
14700 R=M/511.
14800 T=-R
14900 DO 500 L=1,511
15000 T=T+R
15100 KT=T+1
15200 KU=T
15300 TT=T-KU
15400 X1(L)=SLM3(L1)
15500 500 Y1(L)=SLM3(L2)
15600 X1(512)=X(N)
15700 Y1(512)=Y(N)
15800 600 RETURN
15900 END
16000
16100
16200 C GIVEN TWO POINTS (XX1,YY1), (XX2,YY2) AND A SLOPE SS2
16300 C THIS ROUTINE FINDS THE SLOPE SS1
16400 SUBROUTINE S102(XX1,YY1,XX2,YY2,SS1,SS2)
16500 LF=0
16600 IF(SS2.NE.0.) GO TO 5
16700 IF((XX1.LT.XX2.AND.YY1.LT.YY2).OR.
16800 1 (XX2.LT.XX1.AND.YY2.LT.YY1)) SS1=1001.
16900 IF((XX1.LT.XX2.AND.YY1.GT.YY2).OR.
17000 1 (XX2.LT.XX1.AND.YY2.GT.YY1)) SS1=-1001.
17100 GO TO 100
17200 5 T=(YY1-YY2)/SS2
17300 X=XX2+T
17400 Y=YY1
17500 XMAX=XX2
17600 XMIN=XX1
17700 IF(XX2.GE.XX1) GO TO 10
17800 XMAX=XX1
17900 XMIN=XX2
18000 10 IF(X.GT.XMIN.AND.X.LT.XMAX) GO TO 20
18100 T=XX1-XX2
18200 Y=YY2+SS2*T
18300 X=XX1
18400 LF=1
18500 20 IF(LF.EQ.1) GO TO 30
18600 D2=XX1-X
18700 D1=X-XX2
18800 R=D2/D1
18900 XX=XX2
19000 YY=(YY2+R*YY1)/(1.+R)
19100 GO TO 40
19200 30 D2=YY1-Y
19300 D1=Y-YY2
19400 R=D2/D1
19500 YY=YY2
19600 XX=(XX2+R*XX1)/(1.+R)
19700 40 SS1=(YY-YY1)/(XX-XX1)
19800 100 RETURN
19900 END
20000
20100 SUBROUTINE S101
20200 COMMON X(100),Y(100),N,X1(512),Y1(512),S(100),K
20205 DATA SZ/1002001./
20210 SS=S(K)**2
20240 SSS=S(K+1)**2
20270 IF(SS.NE.SZ.AND.SSS.NE.SZ)RETURN
20400 IF(SS.EQ.SZ.AND.X(K).LT.X(K+1).AND.
20500 1 Y(K).GT.Y(K+1)) S(K)=-1001.
20600 IF(S(K)**2.EQ.SZ.AND.X(K).GT.X(K+1).AND.
20700 1 Y(K+1).GT.Y(K)) S(K)=-1001.
20800 IF(SSS.EQ.SZ.AND.X(K).LT.X(K+1).AND.
20900 1 Y(K).GT.Y(K+1)) S(K+1)=-1001.
21000 IF(S(K+1)**2.EQ.SZ.AND.X(K).GT.X(K+1).AND.
21100 1 Y(K).LT.Y(K+1)) S(K+1)=-1001.
21200 100 RETURN
21300 END
21400
22000 FUNCTION SLM2(L)
22100 COMMON X(100),Y(100),N,X1(512),Y1(512),S(100),K
22175 M=K+L
22187 J=M+1
22200 SLM2=(Y(J)-Y(M))/(X(J)-X(M))
22300 END
22400
22600 FUNCTION SLM3(L)
22700 COMMON/SLM/A(2,200),B(2,200),C(2,200),I(200),I1(200),
22750 1 I2(200),KT,TT
22800 SLM3=A(L,KT)*TT**2+B(L,KT)*TT+C(L,KT)
22900 END